Relativistic MHD in dynamical spacetimes: 
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We recently developed a new general relativistic magnetohydrodynamic code with adaptive mesh 
refinement that evolves the electromagnetic (EM) vector potential Ai instead of the magnetic fields 
directly. Evolving Ai enables one to use any interpolation scheme on refinement level boundaries 
and still guarantee that the magnetic field remains divergenceless. As in classical EM, a gauge 
choice must be made when evolving Ai, and we chose a straightforward "algebraic" gauge condi- 
tion to simplify the Ai evolution equation. However, magnetized black hole-neutron star (BHNS) 
simulations in this gauge exhibit unphysical behavior, including the spurious appearance of strong 
magnetic fields on refinement level boundaries. This spurious behavior is exacerbated when matter 
crosses refinement boundaries during tidal disruption of the NS. Applying Kreiss-Oliger dissipation 
to the evolution of the magnetic vector potential Ai slightly weakens this spurious magnetic effect, 
but with undesired consequences. We demonstrate via an eigenvalue analysis and a numerical study 
that zero-speed modes in the algebraic gauge, coupled with the frequency filtering that occurs on 
refinement level boundaries, are responsible for the creation of spurious magnetic fields. We show 
that the EM Lorenz gauge exhibits no zero-speed modes, and as a consequence, spurious magnetic 
effects are quickly propagated away, allowing for long-term, stable magnetized BHNS evolutions. 
Our study demonstrates how the EM gauge degree of freedom can be chosen to one's advantage, 
and that for magnetized BHNS simulations the Lorenz gauge constitutes a major improvement over 
the algebraic gauge. 
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I. INTRODUCTION 

Magnetized fluids in dynamical, strongly curved space- 
times play a central role in many systems of current inter- 
est in relativistic astrophysics. The presence of magnetic 
fields may destroy differential rotation in nascent neu- 
tron stars, form jets and influence disk dynamics around 
black holes, affect collapse of massive rotating stars, etc. 
Many of these systems are promising sources of gravi- 
tational radiation for detection by laser interferometers 
such as LIGO, VIRGO, TAMA, GEO, LCGT, ET, LISA 
and DECIGO. Some may also emit electromagnetic ra- 
diation as gamma-ray bursts, or emission from magne- 
tized disks around black holes in active galactic nuclei 
and quasars. Accurate, self-consistent modeling of these 
systems requires computational schemes capable of si- 
multaneously accounting for magnetic fields, relativistic 
magnetohydrodynamics (MHD), radiation transport and 
relativistic gravitation. 

Over the past several years, we have developed a nu- 
merical scheme in 3+1 dimensions capable of solving 
the coupled system of Einstein's field equations, rel- 
ativistic MHD and Maxwell's equations [l|, 0. Our 
approach is based on the BSSN (Baumgarte-Shapiro- 
Shibata-Nakamura) formalism to evolve the metric [H, H| , 
a high-resolution, shock-capturing (HRSC) scheme to 
handle the fluids, and a constrained-transport (CT) 
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scheme to treat the magnetic induction equation [5| . This 
GRMHD code has been subjected to a rigorous suite of 
numerical tests to check and calibrate its accuracy 0, [|| . 
The code has been applied to explore a number of im- 
portant dynamical scenarios in relativistic astrophysics, 
including the collapse of magnetized, differentially ro- 
tating hypermassive neutron stars to black holes [6|-8|, 
the collapse of rotating stellar cores to neutron stars 9[, 
the collapse of rotating, supermassive stars and massive 
Population III stars to black holes [lf3 |. magnetized bi- 
nary neutron star merger (llj , binary black hole- neutron 
stars (BHNSs) fl2l . fill ] , binary white dwarf- neutron stars 
[Til Ha l , an d the mer ger of binary black holes in gaseous 
clouds [16| and disks |17j . 

Many problems in relativistic astrophysics require nu- 
merical simulations covering a large range of length 
scales. For example, to follow the final merger of a com- 
pact binary system with a total mass M, length scales 
of M/30 or less must be resolved to treat the strong- 
field, near-zone regions reliably. On the other hand, ac- 
curate gravitational wave calculations at length scales 
~ M must be performed far in the weak-field wave-zone 
at radius r > lOOAf. Adaptive mesh refinement (AMR) 
allows for sufficient resolution to be supplied to areas of 
the computational domain as needed, enabling us to re- 
solve both strong- and weak-field domains efficiently. 

One of the most subtle issues in evolving the MHD 
equations is the preservation of the divergenceless con- 
straint (V • B = 0). When evolving the induction equa- 
tion, numerical truncation error leads to violations of 
this constraint, resulting in unphysical plasma trans- 
port orthogonal to the magnetic field, as well as viola- 
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tions of energy and momentum conservation (see e.g., 
Ell)- In simulations using a uniforml y-sp aced grid, 
constrained-transport schemes (see e.g., [H.l20||) are com- 
monly used to maintain the divergenceless constraint. In 
simulations using AMR grids, both CT schemes and the 
hyperbolic divergence-cleaning (HDC) method [HI, [HI 
have been used. However, we find that in the presence of 
black holes, otherwise stable GRMHD evolution schemes 
often become unstable. 

One of the most commonly adopted methods for evolv- 
ing black holes is the moving puncture technique [23J, [24| , 
in which the puncture gauge conditions ensure that the 
spacetime singularity in the black hole interior never ap- 
pears. However, a coordinate singularity is present in 
the computational domain near which accurate numerical 
evolution is difficult to achieve. It has been demonstrated 
that when using the moving puncture gauge, while gauge 
modes can leak out of the BH horizon [25| , any physical 
or constraint violating data in the black hole interior will 
not propagate out of the horizon [261-28] . We find that 
this property is preserved in the presence of hydrody- 
namic matter, but can be challenging to maintain when 
the matter is magnetized. Though we have found that 
the HDC scheme works well with the moving puncture 
technique in the absence of black holes, when BHs are 
present and excision is not applied, we find that, even in 
the Cowling approximation, in which the metric is fixed, 
inaccurate data in the black hole interior can propagate 
out of the horizon and contaminate the solutions. 

For this reason, in developing a moving puncture- 
compatible algorithm for maintaining V • B = 0, we fo- 
cused on CT schemes. This was the approach adopted 
in our earlier, unigrid GRMHD code In AMR simu- 
lations, the conventional CT scheme may be used where 
individual refinement levels do not overlap. However, 
maintaining the divergenceless constraint at refinement 
level boundaries requires that special interpolations be 
performed during prolongation and restriction. Such pro- 
longation/restriction operators have been devised |29l — 
l3lj . but these operators must be fine-tuned to the par- 
ticular AMR implementation, and some have only been 
tested for Newtonian MHD. 

In [2j we developed an alternative, AMR-compatible 
CT scheme. The scheme is based on the CT method 
described in [32|], which recasts the magnetic induction 
equation as an evolution equation for the magnetic vector 
potential. Divergence-free magnetic fields are guaranteed 
as they are computed from the curl of the vector poten- 
tial. The evolution of the vector potential is carried out 
in the same HRSC framework as other hydrodynamic 
variables. This scheme is numerically equivalent to the 
commonly- used, staggered- -B CT scheme of for uni- 
form resolution simulations, but is readily generalized to 
an AMR grid. Unlike the magnetic field, the vector po- 
tential is not constrained, so any interpolation scheme 
may be used during prolongation and restriction, thus 
enabling its use with any AMR algorithm. In addition, 
since the vector potential does not uniquely determine 



the magnetic field, there exists an extra gauge degree of 
freedom that can be used to one's advantage. A straight- 
forward "algebraic" gauge condition was adopted in our 
original formulation, to simplify evolution equation for 
the vector potential (see also [33. l33jV 

In we performed several tests on our new AMR- 
compatible CT scheme. We found that our scheme works 
well even in black-hole spacetimes. Inaccurate data gen- 
erated in the black hole interior stayed inside the horizon. 
Hence our scheme is compatible with the moving punc- 
ture technique. However, when this algebraic gauge con- 
dition is used for magnetized binary black hole-neutron 
star (BHNS) simulations, stable, long-term evolutions are 
not possible. 

In this paper, we exploit another EM gauge condi- 
tion, the Lorenz gauge 3J] , which allows for stable, long- 
term GRMHD evolutions of binary BHNSs. Performing 
eigenvalue analyses, we find that static EM modes are 
present when adopting the algebraic gauge, whereas all 
EM modes are propagating when the Lorenz gauge is 
adopted. As a result, when the algebraic gauge is used, 
its zero-speed modes are inevitably excited and remain on 
the grid. Further, when these modes are interpolated at 
refinement boundaries, inaccuracies in interpolation may 
lead to the appearance of spurious magnetic fields that 
contaminate the solution and lead to a build-up of nu- 
merical error. Given that the Lorenz gauge possesses no 
static modes, any such spurious modes should be prop- 
agated away to the boundary and leave the grid. Thus, 
we expect that the solutions will be improved when using 
the Lorenz gauge. 

Using magnetized BHNS simulations we perform a 
comparison between the Lorenz gauge and the algebraic 
gauge with and without Kreiss-Oliger dissipation (KOD) 
applied on the vector potential evolution. We confirm 
the expected behavior. The BHNS system that forms 
the basis of our comparison is an irrotational binary in 
an initially circular orbit, with a BH:NS mass ratio of 3:1 
and the NS initially seeded with a small poloidal mag- 
netic field. Using these simulations we demonstrate the 
superiority of the Lorenz gauge, and confirm that in the 
algebraic gauge static modes are present and lead to spu- 
rious B-fields which remain on the grid. This spurious 
effect contaminates the solution and is amplified in time, 
leading the simulations to break down shortly after NS 
tidal disruption. In contrast to the algebraic gauge, static 
EM modes are not present in the Lorenz gauge, and any 
spurious fields generated at refinement level boundaries 
are quickly propagated away. Overall we find that the 
Lorenz gauge provides a major improvement for magne- 
tized BHNS simulations. In J35|, we employ the Lorenz 
gauge to perform a detailed study of the effects of mag- 
netic fields in the evolution of binary BHNSs. 

The remaining sections are organized as follows. In 
Sec. |TT] we outline the basic GRMHD equations and pro- 
vide an overview of EM gauge conditions and evolution 
equations. In Sec. IIIII we perform an eigenvalue analysis 
of the EM evolution equations. In Sec. IIV1 we briefly 
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describe the specific implementation of our GRMHD 
scheme. In Sec. |Vj we review the differences between 
our magnetized BHNS simulations in the algebraic and 
Lorenz gauges. We summarize our findings in Sec. IVI1 



II. BASIC EQUATIONS 

The formulation and numerical scheme for our simula- 
tions is basically the same as in our previous work 0, 
[ill ], to which the reader may refer for details. Here we in- 
troduce our notation, summarize our method, and point 
out the latest changes to our numerical technique. Ge- 
ometrized units (G = c = I) are adopted, except where 
stated explicitly. Greek indices denote all four spacetime 
dimensions (0, 1, 2, and 3), and Latin indices imply spa- 
tial parts only (1, 2, and 3). 

We use the 3+1 spacetime decomposition, in which the 
metric has the following form: 



ds 2 



-a 2 dt 2 + ^{dx 1 + ftdt^dx 1 + ft dt) 



where a is the lapse function, /3* the shift vector, and 
7ij the three-metric on spacelike hypersurfaces of con- 
stant time t. We then decompose all evolution equations 
in 3+1 form (see e.g. [36| for a detailed discussion and 
references) . 



A. Gravitational fields 

We adopt the Baumgarte-Shapiro-Shibata-Nakamura 
(BSSN) formalism [3j, |4| in which the metric evolution 
variables are the conformal exponent <fi = ln(7)/12, the 
conformal 3-metric jij — e ^Tij, three auxiliary func- 
tions T l = — 7 y ,j, the trace K of the extrinsic curvature 
Kij, and the trace- free part of the conformal extrinsic 
curvature Aij = e^ 4 ^(Kij — jijK/3). Here, 7 = det(7jj). 
The full spacetime metric <? p „ is related to the three- 
metric 7 MI , by 7 M „ = g^ v + n^n v , where the future- 
directed, timelike unit vector normal to the time 
slice can be written in terms of the lapse a and shift 
f3 l as n M = a _1 (l, —f3 l ). The evolution equations of 
these BSSN variables are given by Eqs. (9)-(13) of [Uj]. 
We adopt the standard puncture gauge conditions: an 
advective "1+log" slicing condition for the lapse and a 
"Gamma- freezing" condition for the shift [37]. The evo- 
lution equations for a and (3 Z are given by Eqs. (2)-(4) 
of [111, with tne V parameter set to w 2.2/M in all BHNS 
simulations presented here, where M is the ADM mass 
of the binary. We add a sixth-order KOD term of the 
form e/6A(Ax 5 d { x 6) +Ay 5 di 6) +Az 5 di 6) )f to reduce high- 
frequency numerical noise associated with AMR refine- 
ment interfaces. Here / stands for all evolved BSSN, 
lapse and shift variables. We choose the strength param- 
eter as e = 0.2. Note that this KOD term converges to 
zero at fifth-order in the grid spacing so that it does not 
alter the convergence properties of our BSSN scheme. 



B. Magnetohydrodynamic fields 

The fundamental MHD variables include the rest-mass 
density po, specific internal energy e, pressure P, four- 
velocity u M , and magnetic field B^ = n l/ F* vfl . Here F*^ v 
is the dual of the Faraday tensor F^ v . Note that is 



purely spatial (B 



L B^/a = 0). We adopt a T-law 



equation of state (EOS) P = (T — l)p e with T = 2, 
which reduces to an n — 1 polytropic law for the initial 
(cold) neutron star matter. The stress-energy tensor is 
given by 



Tp V = (pah + b )u^u^ + I P + — J - b^K , (2) 
where h = 1 + e + P/ po is the specific enthalpy and 



P B v 



An n,,u v 



(3) 



(1) is the magnetic field measured in fluid's comoving frame, 



modulo a factor of 1/v4tt- Here P^ v = g^ u + u^u v and 
b 2 = In the ideal MHD limit, in which the plasma is 

assumed to have perfect conductivity, the Faraday tensor 
can be written as F^ v = \fl~Tr u^^bg. 

In the standard numerical implementation of the MHD 
equations using a conservative scheme, it is useful to in- 
troduce the "conservative" variables p*, Si, f and B l . 
They are defined as 



P* = -VlPon^ , 
Si = - v / 7 7 V n 'V l , 

B l = ^B\ 



(4) 

(5) 
(6) 
(7) 



The evolution equations for p*, Si and f can be derived 
from the conservation of baryon number V Al (p*u /i ) = 
and conservation of energy-momentum V M T Ml/ = 0, giv- 
ing rise to Eqs. (27)-(30) in 0. 



C. Electromagnetic fields 

In the ideal MHD limit, the Maxwell equation 
V y F*^ v — yields the magnetic constraint djB^ = 
and the induction equation dtB 1 + dj(v^ B l — v 1 B' J ) = 0. 
As shown in [2J and [38|, these equations can be rewrit- 
ten by introducing the electromagnetic 4-vector poten- 



tial A u 



A„, with n^Af, 



0. The magnetic 



constraint and induction equations become 



dtAi = e ijk v j B k -di(a$ 



■FA;) 



(8) 
(9) 



where e^ k — n Al e M *-' fe is the 3-dimcnsional Levi-Civita 
tensor. The dynamical variable in this scheme is the 
vector potential. The B-field in Eq. ((9]) and in the MHD 
equations is essentially replaced according to Eq. ©. 
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The EM evolution via the vector potential introduces 
an additional gauge degree of freedom. Therefore, an 
EM gauge choice must be made to close the system of 
equations. 

1. Algebraic gauge 

In 0, 32], the algebraic EM gauge condition 

$ = -/3 j Aj = -n j Aj (10) 
is adopted because it simplifies the evolution of Ai to 

d t A t = <,,!.<■■ ir. (11) 

2. Lorenz gauge 

In this paper, we find that simulations of magnetized 
BHNSs using AMR and the algebraic gauge suffer from 
severe numerical errors and eventually break down due 
to a build-up of errors. 

Improved, stable, long-term BHNS evolutions may be 
achieved by imposing the Lorenz gauge condition 

V^" = 0, (12) 
which yields the evolution equation 

ft(V7*) + fli(«V7^- V7/3 J $) =0 (13) 

for $. 

Our numerical implementation of Eqs. (|8|) and (|9]) 
guarantees numerically identical B l regardless of EM 
gauge in simulations with a uniform-resolution grid. 
However, interpolations performed on Ai at refinement 
boundaries on AMR grids will modify Ai, resulting in 
different B l on and near these boundaries, ultimately 
leading to different stability properties. 

III. ELECTROMAGNETIC GAUGE ANALYSIS 

The differences in the stability properties of the gauges 
considered here can be explained by the fact that the al- 
gebraic gauge possesses a static mode, whereas all modes 
in the Lorenz gauge are propagating. Here we perform 
an eigenvalue analysis of the EM field evolution equa- 
tions with the algebraic and Lorenz gauges and demon- 
strate explicitly the existence of static modes in the for- 
mer gauge. 

A. Evolution Equations 

Substituting Eq. ([5} into Eq. © we find that the evo- 
lution of the vector potential is given by 

d t Ai = - v m d m A l + (v m + p m )diA m - adt® 

-^d i a + A m d i m . ( ' 



Writing Eq. (|T3| as an evolution equation for $ we find 

9 t $ = (3 m d m <S> - ad m A m + j- 1/2 ^d m (y^(3 m ) 
- 1 - 1 ' 2 A m d m {a^) - $d t ln(V7)- 

(15) 

In the limit of short-wavelength, small-amplitude per- 
turbations about a background solution, the EM evolu- 
tion equations decouple from the MHD and BSSN equa- 
tions and can be analyzed separately. While the entire 
coupled system of the BSSN, MHD and Maxwell's equa- 
tions is quasi-linear, the Ai,$ partial differential equa- 
tions alone are linear, as the coefficients that multiply 
the spatial derivatives of Ai and $ are not functions of 
Ai and For our purpose we can apply the hyperbolic- 
ity theorems for linear systems of equations [39|, Ho] , thus 
we need only consider the principal parts of the evolution 
equations and the localization principle. 

B. Algebraic gauge 

In this gauge only the vector potential evolves dynam- 
ically and $ is entirely determined by Ai and the BSSN 
gauge variables. The principal part of the Ai evolution 
equation is given by 

d t Ai ~ -v m d m Ai + v m d l A m . (16) 

Here ~ denotes "the principal part equals" . Replacing 
di by ki, where ki is the spatial component of the wave 
vector, normalized according to j^kikj — 1, we obtain 

d t A t ~ -v m k rn A % + v m hA m = Fi (17) 

We can now construct the characteristic matrix of the 
system: = dFi/dAj. For the system to be strongly 
hyperbolic the eigenvalues of this matrix must be real 
and the matrix must possess a complete set of eigenvec- 
tors [!l|]. A straightforward calculation shows that the 
set of eigenvectors of this system is complete and the 
eigenvalues (or characteristic speeds) are 

Ai = 0, A 2 , 3 = v m k m . (18) 

Hence the system is strongly hyperbolic. However, the 
system possesses a static gauge mode, i.e., a mode that 
does not propagate. In unigrid simulations, this is not an 
issue since our numerical implementation of the Ai evo- 
lution guarantees B l (obtained via Eq. (j8])) to be equiva- 
lent to the standard staggered- B CT scheme [32] ]. which 
evolves B % directly and thus does not suffer from poten- 
tial EM gauge mode problems. In AMR simulations how- 
ever, Ai is interpolated between refinement levels. The 
static modes will remain in the computational domain 
where they were generated, allowing interpolation errors 
between refinement levels to give rise to spurious B l at 
refinement level boundaries, which subsequently propa- 
gate to the interior of the refinement boxes and affect 
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TABLE I. Summary of comparison runs. The third column 
specifies whether KOD is applied, and the fourth column in- 
dicates whether build-up of spurious £?-field is observed. A / 
means "Yes" and a X means "No" . 

Case Name Gauge KOD Spurious B r 

Algebraic Algebraic X / 

Algebraic-KO Algebraic / / 

Lorenz Lorenz X X 



the dynamics of system. In this way, numerical errors 
can build up and lead to unstable behavior. On a related 
note, we point out that the existence of zero modes has 
also been considered as an explanation for the different 
stability properties of different 3+1 formulations of GR 
[43 - |45j . However, the existence of zero speed modes is 
not always associated with unstable behavior. For exam- 
ple the BSSN formulation gives rise to stable numerical 
integrations even though it possesses a zero-speed mode. 

C. Lorenz gauge 

In this gauge both the vector and scalar potentials 
evolve. The principal part of the EM evolution equations 
in the Lorenz gauge is 

8 t A t ~v m d m A z + (v m + Dd t A m - ad*®, 
d t $ ~ f3 m d m $ - ai mn d rn A n . ' ' 

Using the &i — > ki rule and calculating the eigensystem 
of the characteristic matrix we find that the matrix has 
a complete set of eigenvectors and its eigenvalues are 

Ai, 2 = P m k m ± a, X 3 a = v m k m . (20) 

As in the algebraic gauge, all eigenvalues are purely real. 
The system is strongly hyperbolic, only this time all 
modes are non-zero. Thus numerical errors arising from 
refinement boundary interpolations will propagate away. 

This is precisely what we observe in our BHNS sim- 
ulations and in the following sections we show explicit 
numerical results that demonstrate this. 



IV. NUMERICAL METHODS 
A. Initial data 

BHNS quasiequilibrium initial data are constructed as 
described in J46|. The particular BHNS binary chosen 
in our comparison of different EM gauges corresponds 
to Case A of [13]. The system consists of an initially 
nonspinning BH orbiting about an irrotational NS, with 
a 3:1 BHNS mass ratio. The NS has a compaction of 
C = 0.145, and is constructed using an n = 1 (T = 2) 
polytropic EOS. 



Given that interior NS magnetic field strengths and 
configurations are unknown, we choose an initial seed 
poloidal magnetic field via a vector potential of the form 

A t = (-^pS^ + ^p-S» % ) A v (21) 

A v = A b w 2 c max(P - P cut , Q) n » (22) 

where (x c ,y c ,0) is coordinate location of the center of 
mass of the NS, w 2 — (x — x c ) 2 + (y — y c ) 2 , and At, 
n p and P cu t are free parameters. The cutoff pressure 
parameter P cut confines the B-field within the neutron 
star to P > P cu t- The parameter n& determines the de- 
gree of central condensation of the magnetic field. Sim- 
ilar profiles of initial magnetic fields have been used 
in numerical simulations of magnetized accretion disks 
(see, e.g. I47l l48l ) and magnetized compact binaries (see, 
e.g. [H, Hg |5Q|). We set P cut to be 4% of the maxi- 
mum pressure and n&=2. The parameter A b controls the 
strength of the initial magnetic field and can be charac- 
terized by the maximum magnetic field inside the NS, 
which we set here to B = 1.3 x 10 16 G (assuming the NS 
mass is 1.4M Q ). Such seed magnetic fields are too weak 
to significantly perturb the initial quasiequilibrium NS, 
and thereby do not contribute to any significant initial 
gravitational field constraint violations. 

B. Evolution algorithms 

Our methods for evolving the metric and the MHD 
equations are described in [13] and for this reason are 
not presented here. This section focuses solely on our 
treatment of the evolution of the vector potential. 

The magnetic induction equation is evolved via the 
4- vector potential using Eqs. ^ and (|13l) . Ai and B % 
are staggered as in and $ exists on the (i + ,j + ,k + ) 
staggered grid [all other hydrodynamic, BSSN, lapse and 
shift variables are stored at k)], where (i + , j + , k + ) = 
(i + 1/2, j + 1/2, k + 1/2). The term -80^^) in 
Eq. (p~3|) is treated using a second-order upwind scheme, 
and Eq. ([9]) is handled using the finite- volume equations 
similar to Eqs. (63)^(65) in [2j. For example, we evolve 
A x , which is staggered at (i, j + ,k + ), by 

dt{A x ) l , J +,k+ = -£?,j+,k+ ~ ^t( a$ _ ft A i)i+ ,j+,k+ 
-{a*-pA j )*- tj+M ], (23) 

where i~ = i — 1/2, A x is the line-averaged A x defined 
in Eq. (57) of 0, £ l = e i jk B^v h and £ x is the line- 
averaged £ x defined in Eq. (49) of 0- In @ and [H, 
only the first term on the right-hand side of Eq. (f2U| is 
present because the algebraic gauge Eq. pH|) is imposed. 
The term £ x - + k+ is computed using the HRSC scheme 

described in Q and [23]. Values of a, ft and Ai at 
the staggered points (i ± ,j + ,k + ) on the right-hand side 
of Eq. (|23p are computed by simple averaging. Similar 
equations can be derived for A y and A z . 
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TABLE II. Grid configuration used in our simulations. Nah denotes the number of grid points covering the diameter of the 
apparent horizon initially, and Nss denotes the number of grid points covering the smallest diameter of the NS initially. 

Grid Hierarchy (in units of M)' a ' Max. resolution A?ah -Ans 

(203.7, 101.8, 25.46, 12.73, 6.365, 3.183, 1.591) M/31.4 40 82 



There are two sets of nested refinement boxes: one centered on the NS and one on the BH. This column specifies the half 
side length of the refinement boxes centered either on the BH or the NS. 



Ai is staggered to ensure that B fields obtained by 
taking the curl operator on Ai [Eqs. (60)-(62) in Q] are 
numerically equivalent to those obtained from the stan- 
dard, staggered- B constrained transport scheme [20| . It 
is easy to prove that the additional terms in Eq. (l23l) and 
the corresponding terms in A y and A z equations cancel 
exactly after applying the curl operator. The resulting 
numerical values of B l are thus gauge-invariant in un- 
igrid simulations. We have confirmed numerically that 
this is indeed the case. However, in simulations with an 
AMR grid, since interpolations on Ai are performed be- 
tween refinement levels, values of Ai are not the same in 
different EM gauges at the refinement boundaries. The 
resulting B field at the refinement boundaries are also 
different in general but should converge to a unique, true 
solution with increasing resolution in any gauge. 



C. Recovery of primitive variables 

At each timestep, the "primitive variables" po, P, and 
v % are recovered from the "conservative" variables p», f, 
and Si. We perform the inversion by numerical solving 
two nonlinear e quat ions via the Newton- Raphson method 
as described in [51] . using the public code developed by 
Noble et al. [52| , 

Sometimes numerical errors may drive "conservative" 
variables to unphysical values, resulting in unphysical 
primitive variables after inversion (e.g. negative pressure 
or even complex solutions). This usually happens in the 
low-density "atmosphere" or deep inside the BH interior 
where high-accuracy evolutions are difficult to maintain. 
Various techniques have been proposed to handle inver- 
sion failures (see, e.g. [HI). Our approach imposes con- 
straints on the conservative variables to reduce inversion 
failures (for details see [Hj]). 



V. RESULTS 

A. Cases and grid structure 

TableUsummarizes all simulations presented here. The 
only parameters varied in the BHNS simulations shown 
here are the EM gauge condition and the application 
of fourth-order KOD to the Aj evolution of the form 
e/16(Ax 3 di 4) + Ay 3 4 4) + Az 3 di 4) )f (with strength pa- 



rameter set to 0.1). Here / denotes the vector potential. 
Note that this KOD term converges to zero at third-order 
in the grid spacing, so that it does not alter the conver- 
gence properties of our EM evolution scheme. One sim- 
ulation adopts the Lorenz gauge without KOD (named 
"Lorenz" ) , and two simulations adopt the algebraic gauge 
- one with and one without KOD, called "Algebraic-KO" 
and "Algebraic" , respectively. The grid structure used in 
our simulations is presented in Table |TT] and is the same 
for all cases. Our main findings are summarized in Figs.[T] 
andU 



B. Algebraic gauge 

Figured] demonstrates that in the algebraic gauge (left 
two columns), there is an eigenmode of the A^-field which 
does not propagate (as derived in Eq. (fT5)l ) leading to a 
path of Ai fields that trails the orbiting NS. Two nested 
sets of adaptive-mesh refinement boxes surround and co- 
move with the NS and BH, with lower resolution at larger 
distances from these objects (see Table [TTJ) . When refine- 
ment levels cross the Ai trail, the Ai fields are filtered 
to lower or higher frequencies. This filtering leads to 
the appearance of spurious, strong magnetic fields on re- 
finement box boundaries (see Fig. [5]). When the spu- 
rious, strong magnetic fields appear, the magnetic en- 
ergy outside the horizon jumps by about 5% after about 
half an orbit - a strong effect. Kreiss-Oliger dissipation 
helps to damp the spuriously generated magnetic fields, 
resulting in an initial net drop of magnetic energy out- 
side the horizon by 5% after about half an orbit. Cor- 
responding to this drop in magnetic energy, KOD grad- 
ually smooths the initially steep A^-field gradient in the 
NS, ultimately destroying the NS's magnetic field struc- 
ture (inset of middle frame in Fig. [2]) . Just before tidal 
disruption the magnetic energy outside the horizon in the 
Algebraic and the AlgebraicKO cases is about 17% and 
77% higher than the magnetic energy in the Lorenz case, 
respectively. Moreover, when KOD is turned on, the b 2 
contours lag behind the corresponding b 2 contours when 
employing either the algebraic gauge without KOD or 
the Lorenz gauge (compare the insets in Fig. ^ . 
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FIG. 1. Snapshots of |A| (the magnitude of Ai) profiles in the orbital plane of a BHNS binary plotted at selected times 
according to |A| = 10 _3_0 ' 6t , i = 0, . . . , 5, with darker greyscaling for stronger |A|. Notice the Ai trail left behind the orbiting 
NS in the algebraic gauge with or without KOD, and the absence of such trail in the Lorenz gauge. This trail is an immediate 
consequence of the existence of a static mode in the algebraic gauge. When the refinement boxes intersect the Ai trail, spurious 
magnetic fields are generated. The filled black circled near the center denotes the BH's horizon. 
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FIG. 2. Orbital-plane snapshots of the comoving magnetic energy density (fa 2 ) just before NS tidal disruption, b 2 profiles are 
plotted according to b 2 = W~ 4A9 ~ 1A ' /M 2 , i = 0, . . . ,5, with darker greyscaling for larger values. The insets zoom in on the 
NS and demonstrate that the B-field structure is maintained in the outer layers of the NS when using Lorenz gauge, whereas 
it is blown away at this time when the algebraic gauge is used with or without KOD. Note also the spurious B-fields that have 
appeared at refinement level boundaries using the algebraic gauge. These snapshots clearly demonstrate that the algebraic 
gauge with or without KOD generates spurious magnetic fields and destroys the B-field structure of the NS. 



C. Lorenz gauge 

In contrast to the algebraic gauge, there are no zero- 
speed modes in the Lorenz gauge (see Eq. ([20])). and 
no Ai field trail appears behind the NS (third-column 
plots in Fig.[IJ. Immediately after the simulation begins, 
an initial burst of an A4 gauge wave (with a negligibly 
small B field associated with it) emanates from the NS, 
propagates outwards and is partially absorbed by the BH. 
This wave travels at the speed of light, as expected for the 
gauge modes Ai : 2 in Eq. (|20|) . As the NS orbits the BH, 
the Ai field comoves with the plasma. Therefore, strong, 
spurious magnetic fields at refinement boundaries do not 
appear (right frame in Fig. [5]). Also, since KOD is not 
used, the original magnetic field structure of the NS is 
better-maintained (inset of right panel in Fig. [5]) than in 
the algebraic gauge simulations. 



VI. SUMMARY 

In developing an GRMHD algorithm that maintains 
V • B = and is compatible with the moving punc- 
ture technique, we have focused on CT schemes. Such 
schemes guarantee that V • B = is maintained to trun- 
cation error. In Q, we developed an AMR-compatible 
CT scheme, in which the magnetic induction equation 
is recast as an evolution equation for the magnetic vec- 
tor potential. The divergence-free magnetic field is then 
computed via the curl of the vector potential. Unlike the 
magnetic field, the vector potential is not constrained, 
and so any interpolation scheme can be used during pro- 
longation and restriction, thus enabling its use with any 
AMR algorithm. In addition, given that the vector po- 
tential does not uniquely determine the magnetic field, 
there is an extra gauge degree of freedom that can be 



used to one's advantage. 

In we performed several tests on our new AMR- 
compatible CT scheme. Using an algebraic EM gauge 
condition we found that our scheme works well even in 
black-hole spacetimes. Hence our scheme is compatible 
with the moving puncture technique. However, when we 
use this EM gauge condition in the evolution of mag- 
netized binary BHNSs, stable, long-term simulations are 
not possible. 

In this paper, we introduced a new EM gauge con- 
dition, the Lorenz gauge, which allows for stable, long- 
term GRMHD evolutions of binary BHNSs. Performing 
an eigenvalue analysis we found that static EM modes 
are present when adopting the algebraic gauge, whereas 
all modes are propagating when the Lorenz gauge is 
adopted. 

Using magnetized BHNS simulations, we confirmed the 
expected behavior implied by the eigenmode analysis. As 
the magnetized NS orbits the BH, static modes in the al- 
gebraic gauge lead to a trail of nonzero Ai left behind the 
NS. When such static modes from the algebraic gauge are 
interpolated at refinement boundaries, spurious B-ficlds 
are generated which remain on the grid. This spurious 
effect contaminates the solution and is amplified in time, 
forcing the simulations to terminate shortly after the NS 
tidal disruption. In contrast to the algebraic gauge, when 
the Lorenz gauge is adopted, such spurious effects are 
quickly propagated away to the boundary and leave the 
computational domain. Overall we find that the Lorenz 
gauge provides a major improvement over the algebraic 
gauge for mag netized BHNS simulations. In a forthcom- 
ing paper [35| we employ the Lorenz gauge to perform 
a detailed study of the effects of magnetic fields in the 
evolution of binary BHNSs. 

Despite the excellent behavior observed for BHNS sim- 
ulations when using the Lorenz gauge, this gauge choice 
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may not be the optimal one for all situations. For ex- 
ample, we find that the algebraic gauge yields a better 
result in the magnetized Bondi accretion problem with 
radial magnetic fields. This is likely due to the fact that 
the system is stationary and the evolution in the alge- 
braic gauge rigorously preserves the stationary property: 
dtAi = Eijkvi B k = for radial v % and B l . It is therefore 
useful to explore different EM gauge conditions for differ- 
ent problems. Alternatively, one may consider evolving 
B % directly and construct divergence-free interpolation 
schemes between refinement levels [29T - [3lj |. although im- 
plementation of these schemes may be less straightfor- 



ward. 
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